%% ARIMA procedure
% 1. Plot the data. Identify any unusual observations.
% 2. If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.
% 3. If the data are non-stationary: take first differences of the data until the data are stationary.
% 4. Examine the ACF/PACF: Is an AR(p) or MA(q) model appropriate?
% 5. Try your chosen model(s), and use the AICc to search for a better model.
% 6. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
% 7. Once the residuals look like white noise, calculate forecasts.

%%
% load data
close all;
load '../data/lsLNH.mat';

ls1M = val(:,2);
% 1. Plot the data. Identify any unusual observations.
figure(1);
subplot(2,2,1);plot(dates,ls1M);
datetick('x',20); axis tight;
title('Interbank interest rate 1month');

% plot the return of interbank interest rate 
rls1M = ls1M(2:end) - ls1M(1:end-1);
subplot(2,2,2); plot(dates(2:end),rls1M);
datetick('x',20); axis tight;
title('Return of interbank interest rate 1month');


% 2. If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.
[transdat, lambda] = boxcox(ls1M);

subplot(2,2,3); plot(dates(1:end),transdat,'r');
datetick('x',20); axis tight;
title('Boxcox transform interbank interest rate 1month');

rtransdat = transdat(2:end) - transdat(1:end-1);
subplot(2,2,4); plot(dates(2:end),rtransdat,'r');
datetick('x',20); axis tight;
title('return of Boxcox transform interbank interest rate 1month');
%%

% draw the plot to check if residual is normal distribution
figure(2);
subplot(2,2,1);histfit(rls1M,30);title('histogram for interest rate 1M');
subplot(2,2,3);histfit(rtransdat,30);title('histogram for transformed interest rate 1M');
subplot(1,2,2)
% 
% mu=mean(rtransdat);
% sg=std(rtransdat);
x=linspace(-3,3,200);
% x=linspace(mu-4*sg,mu+4*sg,200);
% pdfx=1/sqrt(2*pi)/sg*exp(-(x-mu).^2/(2*sg^2));
pdfx = tpdf(x,1);
plot(x,pdfx);

density
[kurtosis(rls1M) skewness(rls1M) kurtosis(rtransdat) skewness(rtransdat)]

%% draw the t-students
% close all;
% mu=mean(rtransdat);
% sg=std(rtransdat);
% % x=linspace(mu-4*sg,mu+4*sg,200);
% x=linspace(-3,3,400);
% % pdfx=1/sqrt(2*pi)/sg*exp(-(x-mu).^2/(2*sg^2));
% figure(6)
% hold all;
% pdfx = tpdf(x,2);plot(x,pdfx);
% pdfx = tpdf(x,6);plot(x,pdfx);
% pdfx = tpdf(x,12);plot(x,pdfx);
% legend('t=2','t=6','t=12')

% [c,x] = hist(rtransdat,20);
% bar(x,c/trapz(x,c))
close all
ksdensity(rtransdat)

%%
% 3. If the data are non-stationary: take first differences of the data until the data are stationary.

lags = (0:10)';
[h, pvalue, stats] = kpsstest(rtransdat,'lags',lags,'trend',true); %no trend:
%h = 1: unit root nonstationary
%h = 0 -> failed to reject trend stationary -> stationary
result = [lags pvalue stats h]

hY1 = adftest(rtransdat, 'model','ts', 'lags',2)
% hY1 = 1 indicates that there is sufficient evidence to auggest that it is trend stationary
%%
% Use the variance ratio test on al four series to assess whether the series are random walk 
hY1 = vratiotest(rtransdat)

%%
% 4. Examine the ACF/PACF: Is an AR(p) or MA(q) model appropriate?

figure(3);
subplot(2,1,1);autocorr(rtransdat);
subplot(2,1,2);parcorr(rtransdat);

% [h,p,Qstat,crit] = lbqtest(transdat,'Lags',[1, 5,10,15])

%%
% 5. Try your chosen model(s), and use the AICc to search for a better model.

N = size(transdat,1);
bic = zeros(4,4);
d = 1;
nn = 4;
for j = 0: 4
    for i = 0: 4
%         model = arima(j,d,i);
        model = arima('D',1,'ARLags',1:j,'MALags',1:i);
%         model.Distribution = 't';
        [fit,~,logL] = estimate(model,transdat,'print',false);
        [~, bic(j+1,i+1)] = aicbic(logL,i+j+1,N);    
    end         
end


[row, col] = find(bic == min(min(bic)));
bestAR = row - 1;
bestMA = col - 1;
[min(min(bic)) bestAR bestMA]
%%
% 6. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
close all;


model01 = arima('D',1,'ARLags',[],'MALags',[]);
[fit01,EstParamCov01,logL01,info] = estimate(model01,transdat);

model02 = arima('D',1,'ARLags',[2],'MALags',[]);
[fit02,EstParamCov02,logL02,info] = estimate(model02,transdat);

model03 = arima('D',1,'ARLags',[2],'MALags',[2]);
[fit03,EstParamCov03,logL03,info] = estimate(model03,transdat);

model11 = arima('D',1,'ARLags',[2],'MALags',[2],'Variance',garch(0,1));
[fit11,EstParamCov11,logL11,info] = estimate(model11,transdat);

model12 = arima('D',1,'ARLags',[2],'MALags',[],'Variance',garch(0,1));
[fit12,EstParamCov12,logL12,info] = estimate(model12,transdat);

model13 = arima('D',1,'ARLags',[2],'MALags',[],'Variance',garch(1,1));
[fit13,EstParamCov13,logL13,info] = estimate(model13,transdat);

model14 = arima('D',1,'ARLags',[2],'MALags',[2],'Variance',garch(1,1));
[fit14,EstParamCov14,logL14,info] = estimate(model14,transdat);

%%
[aic,bic] = aicbic([logL01, logL01, logL03, logL11,logL12, logL13, logL14],...
                    [1 2 3 4 3 4 5],size(transdat,1))
%%
% model.Distribution = 't';
% [fit,EstParamCov,logL,info] = estimate(model,transdat);

E = infer(fit,transdat);

% figure
% subplot(2,2,1)
% plot(E./sqrt(fit.Variance))
% title('Standardized Residuals')
% subplot(2,2,2)
% qqplot(E)
% subplot(2,2,3)
% autocorr(E)
% subplot(2,2,4)
% parcorr(E)



%%
% 7. Once the residuals look like white noise, calculate forecasts.
% model = model12;
fit = fit03;
nForecast = 50;
[Yf,YMSE] = forecast(fit,nForecast,'Y0',transdat);
UB = Yf + 1.96*sqrt(YMSE);
LB = Yf - 1.96*sqrt(YMSE);


transYf = (lambda.*Yf+1).^(1/lambda)
Yf = transYf;

Nstart = 300;
figure
h1 = plot(transdat(end-(Nstart - 1):end),'Color',[.75,.75,.75]);
hold on
h2 = plot(Nstart:Nstart+nForecast,[transdat(end);Yf],'r','LineWidth',2);
h3 = plot(Nstart:Nstart+nForecast,[transdat(end);UB],'k--','LineWidth',1.5);
plot(Nstart:Nstart+nForecast,[transdat(end);LB],'k--','LineWidth',1.5);
% set(gca,'XTick',1:10:N);
set(gca,'XTickLabel',datestr(dates(1:10:N),17));
legend([h1,h2,h3],'Transformed rate','Forecast',...
       'Forecast Interval','Location','Northwest')
title('Transformed interbank rate Forecast')